SUBROUTINE ATURAN_TRAPESIUM(a, b, n, I, fungsi)
	
	IMPLICIT NONE
	
	DOUBLE PRECISION h, x, sigma, fungsi, a, b, I
	INTEGER r, n
	
	h = (b - a) / n
	x = a
	
	I = fungsi(a) + fungsi(b)
	sigma = 0
	r = 1
	
	do r=1,(n-1)
		x = x + h
		sigma = sigma + 2*fungsi(x)
	end do
	
	I = (I + sigma)*h/2
END

SUBROUTINE ATURAN_TITIKTENGAH(a, b, n, I, fungsi)
	
	IMPLICIT NONE
	
	DOUBLE PRECISION b, a, h, x, sigma, I, fungsi
	INTEGER r, n
	
	h = (b - a) / n
	x = a + h/2
	
	sigma = fungsi(x)
	r = 1
	
	do r=1,(n-1)
		x = x + h
		sigma = sigma + fungsi(x)
	end do
	
	I = sigma * h
	
END

subroutine simpson_1_per_3(a, b, n, I, fungsi)
! n harus genap
	implicit none
	double precision a, b, I, fungsi
	integer*2 n
	double precision h, x, sigma
	integer*2 r
	
	h = (b-a)/n
	x = a
	I = fungsi(a) + fungsi(b)
	sigma = 0
	r = 1
	while (r<=(n-1)) do
		x = x + h
		if (mod(r,2).eq.1) then
			sigma = sigma + 4*fungsi(x)
		else
			sigma = sigma + 2*fungsi(x)
		end if
		r = r + 1
	end do
	I = (I+sigma)*h/3
end

subroutine simpson_3_per_8(a, b, n, I, fungsi)
! n harus kelipatan 3
	implicit none
	double precision a, b, I, fungsi
	integer*2 n
	double precision h, x, sigma
	integer*2 r
	
	h = (b-a)/n
	x = a
	I = fungsi(a) + fungsi(b)
	sigma = 0
	r = 1
	while (r<=(n-1)) do
		x = x + h
		if (mod(r,3).eq.0) then
			sigma = sigma + 2*fungsi(x)
		else
			sigma = sigma + 3*fungsi(x)
		end if
		r = r + 1
	end do
	I = (I+sigma)*h*3/8
end

SUBROUTINE ROMBERG(a,b,n,integral,f)
	IMPLICIT NONE
	DOUBLE PRECISION f,a,b,integral,r(8,8),sum
	INTEGER n,i,j
	EXTERNAL ATURAN_TRAPESIUM
	n = n/2
	if (n.gt.8) then
		n = 8
	end if
	
	r(1,1)=0.5*(b-a)*(f(a)+f(b))
	do i=2,n
		call ATURAN_TRAPESIUM(a,b,2**i-1,r(i,1),f);
	end do
	
	do i = 2,n
		do j = 2,i
			r(i,j) = r(i,j-1)+(r(i,j-1)-r(i-1,j-1))/(4**(j-1)-1)
		end do
	end do
	integral = r(n,n)
END

SUBROUTINE ATURAN_GL(a, b, I, n, fungsi) 

	IMPLICIT NONE
	INTEGER r, n
	DOUBLE PRECISION f(n), fgl, I, fungsi, a, b, t(n), dt
	EXTERNAL fgl
 
	r = 2
	if(n.EQ.2) then
		t(1) = sqrt(3.D0)/3
		t(2) = -sqrt(3.D0)/3
	else 
		if(n.EQ.3) then
			t(1) = 0.774596669
			t(2) = 0
			t(3) = -0.774596669
		else
			if(n.EQ.4) then
				t(1) = sqrt((3.D0 + 2.D0 * sqrt(6.D0 / 5.D0)) / 7.D0)
				t(2) = sqrt((3.D0 - 2.D0 * sqrt(6.D0 / 5.D0)) / 7.D0)
				t(3) = -sqrt((3.D0 - 2.D0 * sqrt(6.D0 / 5.D0)) / 7.D0)
				t(4) = -sqrt((3.D0 + 2.D0 * sqrt(6.D0 / 5.D0)) / 7.D0)
			else
				if(n.EQ.5) then
					t(1) = 0.906179846
					t(2) = 0.538469310
					t(3) = 0
					t(4) = -0.538469310
					t(5) = -0.906179846
				else
					if(n.EQ.6) then
						t(1) = 0.932469514
						t(2) = 0.661209386
						t(3) = 0.238619186
						t(4) = -0.238619186
						t(5) = -0.661209386
						t(6) = -0.932469514
					end if
				end if
			end if
		end if
	end if
	
	do r = 1, n
		f(r) = fgl(fungsi, a, b, t(r))
	end do
	
	if(n.EQ.2) then
		dt = (f(1) + f(2))
	else
		if(n.EQ.3) then	
			f(1) = 0.555555556 * f(1)
			f(2) = 0.888888889 * f(2)
			f(3) = 0.555555556 * f(3)
			dt = f(1) + f(2) + f(3)
		else
			if(n.EQ.4) then
				f(1) = (0.347854845 * f(1))
				f(2) = (0.652145155 * f(2))
				f(3) = (0.652145155 * f(3))
				f(4) = (0.347854845 * f(4))
				dt = f(1) + f(2) + f(3) + f(4)
			else 
				if(n.EQ.5) then
					f(1) = (0.236926885*f(1))
					f(2) = (0.478628670*f(2))
					f(3) = (0.568888889*f(3))
					f(4) = (0.478628670*f(4))
					f(5) = (0.236926885*f(5))
					dt = f(1) + f(2) + f(3) + f(4) + f(5)
				else
					if(n.EQ.6) then
						f(1) = (0.171324492*f(1))
						f(2) = (0.360761573*f(2))
						f(3) = (0.467913935*f(3))
						f(4) = (0.467913935*f(4))
						f(5) = (0.360761573*f(5))
						f(6) = (0.171324492*f(6))	
						dt = f(1) + f(2) + f(3) + f(4) + f(5) + f(6)
					end if
				end if
			end if
		end if		
	end if
	
	I = (b - a) /2 * dt	
END

DOUBLE PRECISION FUNCTION fgl(fungsi, a, b, t) 
	IMPLICIT NONE
	DOUBLE PRECISION fungsi, a, b, t, x
	x = ((a + b) + (b - a) * t) / 2
	fgl = fungsi(x)
END 
